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Abstract 

The paper proposes a control-theoretic framework for verification of numerical software systems, and puts 
forward software verification as an important application of control and systems theory. The idea is to transfer 
Lyapunov functions and the associated computational techniques from control systems analysis and convex op- 
timization to verification of various software safety and performance specifications. These include but are not 
limited to absence of overflow, absence of division-by-zero, termination in finite time, presence of dead-code, and 
certain user-specified assertions. Central to this framework are Lyapunov invariants. These are properly constructed 
functions of the program variables, and satisfy certain properties — resembling those of Lyapunov functions — along 
the execution trace. The search for the invariants can be formulated as a convex optimization problem. If the 
associated optimization problem is feasible, the result is a certificate for the specification. 

Index Terms 

Software Verification, Lyapunov Invariants, Convex Optimization. 

I. INTRODUCTION 

SOFTWARE in safety-critical systems implement complex algorithms and feedback laws that control 
the interaction of physical devices with their environments. Examples of such systems are abundant 
in aerospace, automotive, and medical applications. The range of theoretical and practical issues that arise 
in analysis, design, and implementation of safety-critical software systems is extensive, see, e.g., [26], 
[37] , and [22]. While safety-critical software must satisfy various resource allocation, timing, scheduling, 
and fault tolerance constraints, the foremost requirement is that it must be free of run-time errors. 

A. Overview of Existing Methods 

1) Formal Methods: Formal verification methods are model-based techniques [44], [41], [36] for 

proving or disproving that a mathematical model of a software (or hardware) satisfies a given specification, 
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i.e., a mathematical expression of a desired behavior. The approach adopted in this paper too, falls under 
this category. Herein, we briefly review model checking and abstract interpretation. 

a) Model Checking: In model checking [14] the system is modeled as a finite state transition system 
and the specifications are expressed in some form of logic formulae, e.g., temporal or propositional logic. 
The verification problem then reduces to a graph search, and symbolic algorithms are used to perform 
an exhaustive exploration of all possible states. Model checking has proven to be a powerful technique 
for verification of circuits [13], security and communication protocols [33], [38] and stochastic processes 
[3]. Nevertheless, when the program has non-integer variables, or when the state space is continuous, 
model checking is not directly applicable. In such cases, combinations of various abstraction techniques 
and model checking have been proposed [2], [17], [54]; scalability, however, remains a challenge. 

b) Abstract Interpretation: is a theory for formal approximation of the operational semantics of 
computer programs in a systematic way [15]. Construction of abstract models involves abstraction of 
domains — typically in the form of a combination of sign, interval, polyhedral, and congruence abstractions 
of sets of data — and functions. A system of fixed-point equations is then generated by symbolic for- 
ward/backward executions of the abstract model. An iterative equation solving procedure, e.g., Newton's 
method, is used for solving the nonlinear system of equations, the solution of which results in an inductive 
invariant assertion, which is then used for checking the specifications. In practice, to guarantee finite 
convergence of the iterates, narrowing (outer approximation) operators are used to estimate the solution, 
followed by widening (inner approximation) to improve the estimate [16]. This compromise can be a 
source of conservatism in analysis. Nevertheless, these methods have been used in practice for verification 
of limited properties of embedded software of commercial aircraft [7] . 

Alternative formal methods can be found in the computer science literature mostly under deductive 
verification [32], type inference [45], and dataflow analysis [23]. These methods share extensive similar- 
ities in that a notion of program abstraction and symbolic execution or constraint propagation is present 
in all of them. Further details and discussions of the methodologies can be found in [16], and [41]. 
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2) System Theoretic Methods: While software analysis has been the subject of an extensive body of 
research in computer science, treatment of the topic in the control systems community has been less 
systematic. The relevant results in the systems and control literature can be found in the field of hybrid 
systems [11]. Much of the available techniques for safety verification of hybrid systems are explicitly 
or implicitly based on computation of the reachable sets, either exactly or approximately. These include 
but are not limited to techniques based on quantifier elimination [29], ellipsoidal calculus [27], and 
mathematical programming [5]. Alternative approaches aim at establishing properties of hybrid systems 
through barrier certificates [46], numerical computation of Lyapunov functions [10], [24], or by combined 
use of bi simulation mechanisms and Lyapunov techniques [20], [28], [54], [2]. 

Inspired by the concept of Lyapunov functions in stability analysis of nonlinear dynamical systems 
[25], in this paper we propose Lyapunov invariants for analysis of computer programs. While Lyapunov 
functions and similar concepts have been used in verification of stability or temporal properties of system 
level descriptions of hybrid systems [48], [10], [24], to the best of our knowledge, this paper is the first 
to present a systematic framework based on Lyapunov invariance and convex optimization for verification 
of a broad range of code-level specifications for computer programs. Accordingly, it is in the systematic 
integration of new ideas and some well-known tools within a unified software analysis framework that we 
see the main contribution of our work, and not in carrying through the proofs of the underlying theorems 
and propositions. The introduction and development of such framework provides an opportunity for the 
field of control to systematically address a problem of great practical significance and interest to both 
computer science and engineering communities. The framework can be summarized as follows: 

1) Dynamical system interpretation and modeling (Section II). We introduce generic dynamical system 
representations of programs, along with specific modeling languages which include Mixed-Integer 

Linear Models (MILM), Graph Models, and MIL-over- Graph Hybrid Models (MIL-GHM). 

2) Lyapunov invariants as behavior certificates for computer programs (Section III). Analogous to a 

Lyapunov function, a Lyapunov invariant is a real-valued function of the program variables, and 
satisfies a difference inequality along the trace of the program. It is shown that such functions can 
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be formulated for verification of various specifications. 
3) A computational procedure for finding the Lyapunov invariants (Section IV). The procedure is 

standard and constitutes these steps: (i) Restricting the search space to a linear subspace. (ii) Using 

convex relaxation techniques to formulate the search problem as a convex optimization problem, 

e.g., a linear program [6], semidefinite program [9], [55], or a SOS program [42]. (iii) Using convex 

optimization software for numerical computation of the certificates. 

II. Dynamical System Interpretation and Modeling of Computer Programs 
We interpret computer programs as discrete-time dynamical systems and introduce generic models that 

formalize this interpretation. We then introduce MILMs, Graph Models, and MIL-GHMs as structured 

cases of the generic models. The specific modeling languages are used for computational purposes. 

A. Generic Models 

1 ) Concrete Representation of Computer Programs: We will consider generic models defined by a 
finite state space set X with selected subsets X C X of initial states, and X^ C X of terminal states, 
and by a set-valued state transition function / : X h-> 2 X , such that f(x) C X^^x e X^. We denote 
such dynamical systems by S(X, /,X , 

Definition 1: The dynamical system S(X, f, X , X^) is a C-representation of a computer program 
V, if the set of all sequences that can be generated by V is equal to the set of all sequences X = 
(x(0), x(l), . . . , x(t), . . . ) of elements from X, satisfying 

x(0)ex cx, x(t + i)ef(x(t)) Vt e z+ (l) 

The uncertainty in x(0) allows for dependence of the program on different initial conditions, and the 
uncertainty in / models dependence on parameters, as well as the ability to respond to real-time inputs. 

Example 1: Integer Division (adopted from [44]): The functionality of Program 1 is to compute the 
result of the integer division of dd (dividend) by dr (divisor). A C -representation of the program is 
displayed alongside. Note that if dd > 0, and dr < 0, then the program never exits the "while" loop and 
the value of q keeps increasing, eventually leading to either an overflow or an erroneous answer. The 
program terminates if dd and dr are positive. 



5 



int IntcgcrDivision ( int dd,int dr ) 
{hit q = {0}; int r = {dd}; 
while (r >= dr) 

{ q = q+i; 

r = r — dr; } 
return r; } 



Z = Zn [-32768, 32767] 

L i J 




X = Z 4 




X Q = {(dd, dr, q, r) e X | q = 0, r = dd} 




X x = {(dd,dr,q,r) eX | r < dr} 




( (dd,dr,q+l,r-dr), 


(dd, dr, q, r) e X\X X 


/:(dd,dr,q,r)^| (dd)dw)) 


(dd, dr, q, r) e X x 



Program 1: The Integer Division Program (left) and its Dynamical System Model (right) 
2 ) Abstract Representation of Computer Programs: In a C-representation, the elements of the state 
space X belong to a finite subset of the set of rational numbers that can be represented by a fixed 
number of bits in a specific arithmetic framework, e.g., fixed-point or floating-point arithmetic. When 
the elements of X are non-integers, due to the quantization effects, the set-valued map / often defines 
very complicated dependencies between the elements of X, even for simple programs involving only 
elementary arithmetic operations. An abstract model over-approximates the behavior set in the interest of 
tractability. The drawbacks are conservatism of the analysis and (potentially) undecidability. Nevertheless, 
abstractions in the form of formal over- approximations make it possible to formulate computationally 
tractable, sufficient conditions for a verification problem that would otherwise be intractable. 

Definition 2: Given a program V and its C-representation S (X, /, X , X^), we say that S (X, /, X , X^ 
is an ^-representation, i.e., an abstraction of V, if X C X, X C X , and f(x) C f(x) for all x E X, 
and the following condition holds: 

Iconic^. (2) 



Thus, every trajectory of the actual program is also a trajectory of the abstract model. The definition 
of Xoo is slightly more subtle. For proving Finite-Time Termination (FTT), we need to be able to infer 
that if all the trajectories of S eventually enter X^, then all trajectories of S will eventually enter X^. 
It is tempting to require that X^ C X^, however, this may not be possible as X^ is often a discrete set, 
while Xqo is dense in the domain of real numbers. The definition of X^ as in (2) resolves this issue. 
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Construction of S(X, /, X , X^) from S(X, /, X , X^) involves abstraction of each of the elements 
X, /, X , Xqo in a way that is consistent with Definition 2. Abstraction of the state space X often involves 
replacing the domain of floats or integers or a combination of these by the domain of real numbers. 
Abstraction of X or X^ often involves a combination of domain abstractions and abstraction of functions 
that define these sets. Semialgebraic set-valued abstractions of some commonly-used nonlinearities are 
presented in Appendix I. Interested readers may refer to [49] for more examples including abstractions 
of fixed-point and floating point operations. 
B. Specific Models of Computer Programs 

Specific modeling languages are particularly useful for automating the proof process in a computational 
framework. Here, three specific modeling languages are proposed: Mixed-Integer Linear Models (MILM), 
Graph Models, and Mixed-Integer Linear over Graph Hybrid Models (MIL-GHM). 

I) Mixed-Integer Linear Model (MILM): Proposing MILMs for software modeling and analysis is 

motivated by the observation that by imposing linear equality constraints on boolean and continuous 

variables over a quasi-hypercube, one can obtain a relatively compact representation of arbitrary piecewise 

affine functions defined over compact polytopic subsets of Euclidean spaces (Proposition 1). The earliest 

reference to the statement of universality of MILMs appears to be [39], in which a constructive proof is 

given for the one-dimensional case. A constructive proof for the general case is given in [49]. 
Proposition 1: Universality of Mixed-Integer Linear Models. Let / : X n- R" be a piecewise affine 

map with a closed graph, defined on a compact state space X C [—1, l] n , consisting of a finite union of 

compact poly topes. That is: 

/ (x) G 2AiX + 2Bi subject to x e X u i (1, N) 

where, each Xi is a compact polytopic set. Then, / can be specified precisely, by imposing linear 
equality constraints on a finite number of binary and continuous variables ranging over compact intervals. 
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Specifically, there exist matrices F and H, such that the following two sets are equal: 

G 1 = {(x,f(x)) \xeX} 

G 2 = {(x,y) \F[ X w v if =y, H[ x w v if = 0, (w, v) e [-1, 1]" 1 " x {-1, 1}""} 

Mixed Logical Dynamical Systems (MLDS) with similar structure were considered in [4] for analysis 
of a class of hybrid systems. The main contribution here is in the application of the model to software 
analysis. A MIL model of a computer program is defined via the following elements: 

1) The state space X C [—1, l] n . 

2) Letting n e = n + n w + n v + 1, the state transition function / : X h-> 2 X is defined by two matrices 
F, and if of dimensions n-by-n e and nn-by-n e respectively, according to: 

f(x) e w v 1] T I H[ XWV 1] T = 0, (w,v) e [-l,lp x {-l,l} n "|. (3) 

3) The set of initial conditions is defined via either of the following: 

a) If X is finite with a small cardinality, then it can be conveniently specified by its elements. 

We will see in Section IV that per each element of X , one constraint needs to be included 
in the set of constraints of the optimization problem associated with the verification task. 

b) If X is not finite, or \X \ is too large, an abstraction of X can be specified by a matrix 
H e R n H xn e which defines a union of compact poly topes in the following way: 

X = {x E X \ H [x w v if = 0, (w,v) e [-1,1]"" x{-l,l}"*}. (4) 

4) The set of terminal states X^ is defined by 

x OQ = {xeX \ H[ x w v if ^0, \/we [-i,ip, W e {-l, i} n "}. (5) 

Therefore, S(X, /,X ,-^oo) is wei l defined. A compact description of a MILM of a program is either 

of the form S (F, H, H ,n,n w ,n v ) , or of the form S (F, H, X ,n,n w ,n v ). The MILMs can represent a 

broad range of computer programs of interest in control applications, including but not limited to control 

programs of gain scheduled linear systems in embedded applications. In addition, generalization of the 

model to programs with piecewise affine dynamics subject to quadratic constraints is straightforward. 
Example 2: A MILM of an abstraction of the IntegerDivision program (Program 1: Section II- A), with 
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all the integer variables replaced with real variables, is given by S (F, H, H , 4, 3, 0) , where 



H 



H 



1 




-2 





1 

-2 





-1 



10 1 

1 1 



2 
-2 




-21001 
10 1 
1 1 



1 

1 

1 1/M 

-1 1 



Here, M is a scaling parameter used for bringing all the variables within the interval [—1,1]. 
2) Graph Model: Practical considerations such as universality and strong resemblance to the natural 

flow of computer code render graph models an attractive and convenient model for software analysis. 

Before we proceed, for convenience, we introduce the following notation: P r . (i, x) denotes the projection 

operator defined as P r (i, x) = x, for all % G ZU {n} , and all x G M. n . 

A graph model is defined on a directed graph G (Af, S) with the following elements: 

1) A set of nodes Af = {0} U {1, . . . , m} U {x} . These can be thought of as line numbers or code 

locations. Nodes and x are starting and terminal nodes, respectively. The only possible transition 

from node n is the identity transition to node tx . 

2) A set of edges £ = k) \ % G Af, j G O (i)} , where the outgoing set O (i) is the set of all 

nodes to which transition from node i is possible in one step. Definition of the incoming set X {%) 
is analogous. The third element in the triplet (i, j, k) is the index for the A;th edge between i and 

j, and An = {k | (i,j,k) G £} . 

3) A set of program variables x t G C R, I g Z (1, n) . Given Af and n, the state space of a graph 

model is X = Af x Q n . The state x — (i,x) of a graph model has therefore, two components: The 

discrete component % G Af, and the continuous component x G £l n C R n . 

4) A set of transition labels T ^ assigned to every edge (i, j, k) G 8, where T ^ maps x to the set T ^x = 



{7j. (x,w,v) | (x,w,v) G S%}, where (w,v) G [-1, lf w x {-1, l} n " , and T% : R n+n * 



+n v 



is a polynomial function and is a semialgebraic set. If T jfi is a deterministic map, we drop 



and define T 



T*(x). 



5) A set of passport labels U.^ assigned to all edges (i, j, k) G £, where U.^ is a semialgebraic set. A 
state transition along edge (i, j, k) is possible if and only if x G n^. 
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6) A set of semialgebraic invariant sets X { C f2 n , i e Af are assigned to every node on the graph, 

such that P r (i,x) G Xj. Equivalently, a state x — (i,x) satisfying x G X\Xi is unreachable. 
Therefore, a graph model is a well-defined specific case of the generic model S(X, f, Xo,Xoo), with 
X = M x IT, X = {0} x X 9 , = {n} x X K and / : X h-> 2 X defined as: 

/(x) = /(i,a;) = {0',r* i a;) | jeO(i), le^ni,}. (6) 
Conceptually similar models have been reported in [44] for software verification, and in [1], [12] 
for modeling and verification of hybrid systems. Interested readers may consult [49] for further details 
regarding treatment of graph models with time-varying state-dependent transitions labels which arise in 
modeling operations with arrays. 
Remarks 

- The invariant set of node contains all the available information about the initial conditions of the 

program variables: P r (0,x) G X$. 

- Multiple edges between nodes enable modeling of logical "or" or "xor" type conditional transitions. 

This allows for modeling systems with nondeterministic discrete transitions. 

k 

- The transition label T ^ may represent a simple update rule which depends on the real-time input. 
For instance, if T = Ax + Bw, and S = R n x [-1, 1] , then x H> {Ax + Bw | w G [-1, 1]} . 
In other cases, T j{ may represent an abstraction of a nonlinearity. For instance, the assignment 
x i — y sin (x) can be abstracted by x >->■ {T (x, w) \ (x, w) G S} (see Eqn. (46) in Appendix I). 

Before we proceed, we introduce the following notation: Given a semialgebraic set II, and a polynomial 

function r : W 1 H> R n , we denote by II (r) , the set: n(r) = {x | r(i) G n}. 

a) Construction of Simple Invariant Sets: Simple invariant sets can be included in the model if they 

are readily available or easily computable. Even trivial invariants can simplify the analysis and improve 

the chances of finding stronger invariants via numerical optimization. 

- Simple invariant sets may be provided by the programmer. These can be trivial sets representing 
simple algebraic relations between variables, or they can be more complicated relationships that 
reflect the programmer's knowledge about the functionality and behavior of the program. 
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- Invariant Propagation: Assuming that Tf- are deterministic and invertible, the set 

is an invariant set for node %. Furthermore, if the invariant sets Xj are strict subsets of f2 n for all 
j Gl (i) , then (7) can be improved. Specifically, the set 

U 4((^)- 1 )nx j ((7;5)- 1 ) (8) 

is an invariant set for node %. Note that it is sufficient that the restriction of T% to the lower 
dimensional spaces in the domains of U.^ and Xj be invertible. 

- Preserving Equality Constraints: Simple assignments of the form : x t (-)■ / (y rn ) result in invariant 
sets of the form Xi = {x \ xi — f (y m ) = 0} at node i, provided that does not simultaneously 
update y rn . Formally, let Tf- be such that (T^x)i—xi is non-zero for at most one element / G Z (1, n) , 
and that {T^x)i is independent of Xf. Then, the following set is an invariant set at node i : 

*i= U {x\ [T*-/]x = 0} 

jei(i), keAij 

3) Mixed-Integer Linear over Graph Hybrid Model (MIL-GHM): The MIL-GHMs are graph models 
in which the effects of several lines and/or functions of code are compactly represented via a MILM. As 
a result, the graphs in such models have edges (possibly self-edges) that are labeled with matrices F and 
H corresponding to a MILM as the transition and passport labels. Such models combine the flexibility 
provided by graph models and the compactness of MILMs. An example is presented in Section V. 
C. Specifications 

The specification that can be verified in our framework can generically be described as unreachability 

and finite-time termination. 

Definition 3: A Program V = S(X, f, X , X^) is said to satisfy the unreachability property with 

respect to a subset X_ C X, if for every trajectory X = x (■) of (1), and every t e Z + , x(t) does 

not belong to X_. A program V = S(X, /, X^X^) is said to terminate in finite time if every solution 

A" = x (•) of (1) satisfies x(t) E Xoo for some t E Z + . 

Several critical specifications associated with runtime errors are special cases of unreachability. 
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1 ) Overflow: Absence of overflow can be characterized as a special case of unreachability by defining: 

X_ = [x E X \ ||« > 1) ol = diagjaj}} 

where on > is the overflow limit for variable i. 

2) Out- of -Bounds Array Indexing: An out-of-bounds array indexing error occurs when a variable 

exceeding the length of an array, references an element of the array. Assuming that x\ is the corresponding 
integer index and L is the array length, one must verify that x t does not exceed L at location i, where 
referencing occurs. This can be accomplished by defining X_ = {(i,x) E X | \x t \ > L} over a graph 

model and proving that X_ is unreachable. This is also similar to "assertion checking" defined next. 

3) Program Assertions: An assertion is a mathematical expression whose validity at a specific location 

in the code must be verified. It usually indicates the programmer's expectation from the behavior of the 
program. We consider assertions that are in the form of semialgebraic set memberships. Using graph 
models, this is done as follows: 

at location % : assert x E Ai =>- define X_ = {(«, x) E X \ x E X\Ai} , 
at location % : assert x ^ A { =>- define X_ = {(i, x) E X \ x E A { } . 
In particular, safety assertions for division-by-zero or taking the square root (or logarithm) of positive 

variables are standard and must be automatically included in numerical programs (cf. Sec. III-A, Table I). 

4) Program Invariants: A program invariant is a property that holds throughout the execution of the 

program. The property indicates that the variables reside in a semialgebraic subset Xj C X. Essentially, 
any method that is used for verifying unreachability of a subset X_ C X, can be applied for verifying 
invariance of X/ by defining X_ = X\Xj, and vice versa. 
D. Implications of the Abstractions 

For mathematical correctness, we must show that if an ^4-representation of a program satisfies the 

unreachability and FTT specifications, then so does the C-representation, i.e., the actual program. This is 

established in the following proposition. The proof is omitted for brevity but can be found in [49]. 
Proposition 2: Let S(X, /, X ,Xoo) be an ^4-representation of program V with C-representation 

S(X, f, X ,Xoo). Let X_ C X and X_ C X be such that X_ C X_. Assume that the unreachability 
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property w.r.t. X_ has been verified for S. Then, V satisfies the unreachability property w.r.t. X_. 
Moreover, if the FTT property holds for S, then V terminates in finite time. 

Since we are not concerned with undecidability issues, and in light of Proposition 2, we will not 
differentiate between abstract or concrete representations in the remainder of this paper. 

III. Lyapunov Invariants as Behavior Certificates 

Analogous to a Lyapunov function, a Lyapunov invariant is a real-valued function of the program 

variables satisfying a difference inequality along the execution trace. 

Definition 4: A (9, fi)- Lyapunov invariant for S(X, f, X , X^) is a function V:1h1 such that 

V(x+)-6V(x)<-n Vx G X, x + G f (x) : x £ X w . (9) 

where {9, if) G [0, oo) x [0, oo). Thus, a Lyapunov invariant satisfies the difference inequality (9) along 
the trajectories of S until they reach a terminal state X^. 

It follows from Definition 4 that a Lyapunov invariant is not necessarily nonnegative, or bounded from 
below, and in general it need not be monotonically decreasing. While the zero level set of V defines an 
invariant set in the sense that V (x k ) < implies V (x k+l ) < 0, for all / > 0, monotonicity depends on 
9 and the initial condition. For instance, if V (x ) < 0, Vx G X , then (9) implies that V (x) < along 
the trajectories of <S, however, V (x) may not be monotonic if 9 < 1, though it will be monotonic for 
9 > 1. Furthermore, the level sets of a Lyapunov invariant need not be bounded closed curves. 

Proposition 3 (to follow) formalizes the interpretation of Definition 4 for the specific modeling lan- 
guages. Natural Lyapunov invariants for graph models are functions of the form 

V(x) = V (i, x) = Oi (x) , i G N, (10) 

which assign a polynomial Lyapunov function to every node % G M on the graph G (A/", £) . 

Proposition 3: Let S (F,H, X ) and properly labeled graph G (J\f, 8) be the MIL and graph 

models for a computer program V. The function V : [—1, 1]™ >->■ R is a (9, /i)-Lyapunov invariant for V 

if it satisfies: 

V(Fx e ) - 9V (x) < -pi, V x e ) G [-1, If x S, 
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where 

E = {(x,w,v,l) | H[ x w v if = 0, (w,v) E [-1,1]"™ x {-1,1}""}. 
The function V : A/"xIR n h-> R, satisfying (10) is a (0, yu)-Lyapunov invariant for P if 

- ^ (x) < -n, V(i,j,k)e£, (x,x + )e(x i nn k ji )xT k ji x. (11) 

Note that a generalization of (9) allows for and fx to depend on the state x, although simultaneous 
search for 9 (x) and V (x) leads to non-convex conditions, unless the dependence of 9 on x is fixed 
a-priori. We allow for dependence of 9 on the discrete component of the state in the following way: 

a 3 {x + )-9%a l {x)<-^ V(i,j,k)e£, (x,x+) e pQ n n^) x T^x (12) 
A. Behavior Certificates 

1) Finite-Time Termination (FTT) Certificates: The following proposition is applicable to FTT analysis 
of both finite and infinite state models. 

Proposition 4: Finite-Time Termination. Consider a program P, and its dynamical system model 

S{X, /, X ,Xoo). If there exists a (9, /x)-Lyapunov invariant V : X i->- R, uniformly bounded on X\Xoo, 

satisfying (9) and the following conditions 

V (x) < -77 < 0, VxG X (13) 
/i + (fl-l)||V|| oo >0 (14) 
max (/i, 77) > (15) 
where || = sup V (x) < 00, then P terminates in finite time, and an upper-bound on the number 
of iterations is given by 

f log + {9 -l)||V||J- log (/x) 



log0 



, 0^1, /i>0 



log (IMP -log fo) a „ n (16) 

> ,0 = 1 



Proof: The proof is presented in Appendix EL 
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When the state-space X is finite, or when the Lyapunov invariant V is only a function of a subset of 
the variables that assume values in a finite set, e.g., integer counters, it follows from Proposition 4 that V 
being a (9, /i)-Lyapunov invariant for any 9 > 1 and \x > is sufficient for certifying FTT, and uniform 
boundedness of V need not be established a-priori. 

Example 3: Consider the integerDivision program presented in Example 1. The function V : X \-> R, 
defined according to V : (dd, dr, q, r) h-> r is a (1, dr) -Lyapunov invariant for IntegerDivision: at every 
step, V decreases by dr > 0. Since X is finite, the program IntegerDivision terminates in finite time. This, 
however, only proves absence of infinite loops. The program could terminate with an overflow. 

2) Separating Manifolds and Certificates of Boundedness: Let V be a Lyapunov invariant satisfying 

def 

(9) with 9 = 1. The level sets of V, defined by C r {V) = {x e X : V(x) < r}, are invariant with respect 
to (1) in the sense that x(t + 1) e C r (V) whenever x(t) e However, for r = 0, the level sets 

£ r (V) remain invariant with respect to (1) for any nonnegative 9. This is an important property with 
the implication that 9 = 1 (i.e., monotonicity) is not necessary for establishing a separating manifold 
between the reachable set and the unsafe regions of the state space (cf. Theorem 1). 

Theorem 1: Lyapunov Invariants as Separating Manifolds. Let V denote the set of all (0, //)- 
Lyapunov invariants satisfying (9) for program V = S(X, f, X , X^). Let I be the identity map, and for 
he {f,I} define 

h- 1 (X_) = {x G X\h (x) n X_ ^ 0} . 
A subset X_ C X, where X_ n X = can never be reached along the trajectories of V, if there exists 
V E V satisfying 

supV(x) < inf V(x) (17) 
and either 9 = 1, or one of the following two conditions hold: 

(I) 9 < 1 and inf V(x) > 0. (18) 

xeh-i-(X-) 

(II) 6" > 1 and supV{x)<0. (19) 

xex 



15 



Proof: The proof is presented in Appendix EL ■ 
The following corollary is based on Theorem 1 and Proposition 4 and presents computationally 

implementable criteria for simultaneously establishing FTT and absence of overflow. 

Corollary 1: Overflow and FTT Analysis Consider a program V, and its dynamical system model 

S(X, /, Xo^Xoo). Let a > be a diagonal matrix specifying the overflow limit, and let X_ = {x G 

x I ll a ~ 1;r lloo > !}■ Let <? e NU{oo} ,he{f,I}, and let the function V : X ^ R be a (0, /i)-Lyapunov 

invariant for 5 satisfying 

V (x) < Vie X . (20) 

F (x) > supjlla- 1 /!^)^ - l} ViGl. (21) 

Then, an overflow runtime error will not occur during any execution of V. In addition, if // > and 

/i + > 1, then, P terminates in at most T u iterations where T u = /i -1 if 9 = 1, and for 7^ 1 we have: 

T _ log (n+{9-l)\\V\\J-\og^ l g(^ + g-l)-l g/i 

log^ S logtf K } 

where II V = sup | V (a;) | . 
xexxix-uxov} 

Proof: The proof is presented in Appendix EL ■ 
Application of Corollary 1 with h — f typically leads to much less conservative results compared with 

h — I, though the computational costs are also higher. See [49] for remarks on variations of Corollary 1 

to trade off conservativeness and computational complexity. 

a) General Unreachability and FTT Analysis over Graph Models: The results presented so far in 

this section (Theorem 1, Corollary 1, and Proposition 4) are readily applicable to MILMs. These results 

will be applied in Section IV to formulate the verification problem as a convex optimization problem. 

Herein, we present an adaptation of these results to analysis of graph models. 

Definition 5: A cycle C m on a graph G (Af, S) is an ordered list of m triplets (ni, n 2 , k±) , (n 2 , n 3 , k 2 ) , 

(n m , n m+ i, k m ) , where n m+1 = n 1: and (nj : rij +1: kj) G S, Vj G Z(l,m). A simple cycle is a cycle 

with no strict sub-cycles. 
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Corollary 2: Unreachability and FTT Analysis of Graph Models. Consider a program V and its 
graph model G (Af, S) . Let V (i, x) = Oi (x) be a Lyapunov invariant for G (Af, E) , satisfying (12) and 

cr (x) < 0, Wxe X 9 (23) 

and either of the following two conditions: 

(I) : a t (x) > 0, Vx G X t n X,_, i G A/"\ {0} (24) 
(II) : a, (x) > 0, Wxe X 3 n T _1 , i G jV\ {0} , j G I (i) , T G {T* } (25) 

where 

T- 1 = {xe Xi\T (x) n ^ 0} 

Then, "P satisfies the unreachability property w.r.t. the collection of sets % G Af\ {0} . In addition, 
if for every simple cycle C G G, we have: 

(0(C)-l)||<7(C)|| oo + Ai(C)>O, and/i(C)>0, and ||<t (C)^ < oo, (26) 

where 

0(C) = II ^, /"(C) = max 4, |k(C)|| 0O = max sup M^l (27) 
then V terminates in at most T u iterations where 

log ((0 (C) - 1) ||<t (C) IU + A* (C)) - log a* (C) , ^ ||(7 (C)|L 



^ " ~ log 0(C) " ^ //(C) 

ceG:e(c)^i 6 v ; ceG:e(C)=i ^ v ; 



Proof: The proof is presented in Appendix II. ■ 
For verification against an overflow violation specified by a diagonal matrix a > 0, Corollary 2 is 
applied with X_ = {x G R n | Ha -1 ^!!^ > 1}- Hence, (24) becomes <7j (x) > p (x) (Wa^xW — 1), 
\/x G Xi, i G A/"\ {0} , where p (x) > 0. User-specified assertions, as well as many other standard safety 
specifications such as absence of division-by-zero can be verified using Corollary 2 (See Table I). 

- Identification of Dead Code: Suppose that we wish to verify that a discrete location % G Af\ {0} in 
a graph model G (Af, £) is unreachable. If a function satisfying the criteria of Corollary 2 with Xj_ = W 1 
can be found, then location % can never be reached. Condition (24) then becomes cr, (x) > 0, Vi 6 M". 
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TABLE I 

Application of Corollary 2 to the verification of various safety specifications. 





apply Corollary 2 with: 


At location i'. 


assert x € X a 


=^ 


:= {x e R n | x € K n \X Q } 


At location i: 


assert x ^ X 


=> 


Xi_ := {i e K" | i G X a } 


At location i: 


(expr.)/cc 


=> 


Xj_ := {i £ R" | i„ = 0} 


At location i: 


2 \/ Xo 


=> 


:= {i€M n |i < 0} 


At location i: 


log (x ) 


=> 


Xi_ := {x e R n \ x o <0} 


At location i: 


dead code 


=> 


:= R n 



Example 4: Consider the following program 



void ComputeTurnRate (void) 

L0 : {double x = {0}; double y = {*PtrToY}; 

LI : while (1) 

L2 : { y = *PtrToY; 

L3 : x= (5 * sin(y) + l)/3; 

L4: ifx>-l{ 

L5 : x = x+ 1.0472; 

L6 : TurnRate = y/x; } 

L7 : else { 

L8 : TurnRate = 100 * y/3.1416 }} 




Program 3. 



Graph of an abstraction of Program 3 



Note that x can be zero right after the assignment x = (5sin(y) + l)/3. However, at location L6, x 
cannot be zero and division-by-zero will not occur. The graph model of an abstraction of Program 3 
is shown next to the program and is defined by the following elements: T 65 : x M- x + 1.0472, and 
T 4 i : x i-> [—4/3, 2] . The rest of the transition labels are identity. The only non-universal passport labels 
are n 54 and n 84 as shown in the figure. Define 

cr 6 (x) = -x 2 - lOOx +1, <x 5 (x) = -(x + 1309/1250) 2 - lOOx - 2543/25 
00 (x) = Oi (x) = cr 4 (x) = cr 8 (x) = —x 2 + 2x — 3. 



It can be verified that V (x) = a { (x) is a (8, l)-Lyapunov invariant for Program 3 with variable rates: 
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# 65 = 1, and % = OV(ij') ^ (6,5). Since 

—2 = sup (T ( x ) < inf °6 ( x ) — 1 

xex xeX - 

the state (6,x = 0) cannot be reached. Hence, a division by zero will never occur. We will show in the 
next section how to find such functions in general. 

IV. Computation of Lyapunov Invariants 
It is well known that the main difficulty in using Lyapunov functions in system analysis is finding 
them. Naturally, using Lyapunov invariants in software analysis inherits the same difficulties. However, 
the recent advances in hardware and software technology, e.g., semi-definite programming [18], [53], 
and linear programming software [19] present an opportunity for new approaches to software verification 
based on numerical optimization. 
A. Preliminaries 

1) Convex Parameterization of Lyapunov Invariants: The chances of finding a Lyapunov invariant are 
increased when (9) is only required on a subset of X\Xoo. For instance, for 9 < 1, it is tempting to 
replace (9) with 

V (x+) - 9V (x) < -/i, Vx G XVXoo :V(x)<l, x + e / (a;) (28) 

In this formulation V is not required to satisfy (9) for those states which cannot be reached from X . 
However, the set of all functions satisfying (28) is not convex and finding a solution for 

(28) is typically much harder than (9). Such non-convex formulations are not considered in this paper. 

The first step in the search for a function V : X \-> R satisfying (9) is selecting a finite-dimensional 
linear parameterization of a candidate function V: 

n 

V (x) = V T (x) = r k V k (x), r = (t*)£ =1 , r k e R, (29) 

k=l 

where 14 : X \-t R are fixed basis functions. Next, for every r = (T k )^ =1 let 



<f>( T ) = v , v maX „ , V t(x+) - 0V T {x), 

xex\x x , x+ef(x) 

(assuming for simplicity that the maximum does exist). Since (■) is a maximum of a family of linear 
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functions, (•) is a convex function. If minimizing (■) over the unit disk yields a negative minimum, 
the optimal r* defines a valid Lyapunov invariant V T *(x). Otherwise, no linear combination (29) yields 
a valid solution for (9). 

The success and efficiency of the proposed approach depend on computability of (•) and its subgra- 
dients. While (■) is convex, the same does not necessarily hold for V T (x+) — 9V T (x). In fact, if X\Xoo 
is non-convex, which is often the case even for very simple programs, computation of (•) becomes a 
non-convex optimization problem even if V T (x+) — V T (x) is a nice (e.g. linear or concave and smooth) 
function of x. To get around this hurdle, we propose using convex relaxation techniques which essentially 
lead to computation of a convex upper bound for (t). 

2) Convex Relaxation Techniques: Such techniques constitute a broad class of techniques for construct- 
ing finite-dimensional, convex approximations for difficult non-convex optimization problems. Some of 
the results most relevant to the software verification framework presented in this paper can be found 
in [31] for SDP relaxation of binary integer programs, [34] and [40] for SDP relaxation of quadratic 
programs, [56] for 5-Procedure in robustness analysis, and [43], [42] for sum-of-squares relaxation in 
polynomial non-negativity verification. We provide a brief overview of the latter two techniques. 

a) The S -Procedure : The 5-Procedure is commonly used for construction of Lyapunov functions 
for nonlinear dynamical systems. Let functions 0j : X i-> R, % e Z (0, m) , and ipj : X i->- R, j e Z (1, n) 
be given, and suppose that we are concerned with evaluating the following assertions: 

(I): 0o (x) > 0, \/x e {x e X \ (pi (x) > 0, ipj (x)=0, ieZ (1, m) , j e Z (1, n)} (30) 

m n 

(II): 3ri e R + , 3fij e R, such that O (x) > (x) + fij^j (x) . (31) 

i=\ j=i 

The implication (II) — > (I) is trivial. The process of replacing assertion (I) by its relaxed version (II) is 
called the 5-Procedure. Note that condition (II) is convex in decision variables and fj,j. The implication 
(I) — > (II) is generally not true and the 5-Procedure is called lossless for special cases where (I) and (II) 
are equivalent. A well-known such case is when m — 1, n — 0, and O , 0i are quadratic functionals. A 
comprehensive discussion of the 5-Procedure as well as available results on its losslessness can be found 



20 

in [21]. Other variations of 5-Procedure with non-strict inequalities exist as well. 

b) Sum-of-Squares (SOS) Relaxation : The SOS relaxation technique can be interpreted as the 
generalized version of the 5-Procedure and is concerned with verification of the following assertion: 

fj (x) > 0, Vj G J, g k (x) ^0,\/ke K, h x (x) =0,VI6L^ -/„ (x) > 0, (32) 

where /_,-, g k , hi are polynomial functions. It is easy to see that the problem is equivalent to verification 
of emptiness of a semialgebraic set, a necessary and sufficient condition for which is given by the Posi- 
tivstellensatz Theorem [8]. In practice, sufficient conditions in the form of nonnegativity of polynomials 
are formulated. The non-negativity conditions are in turn relaxed to SOS conditions. Let E [y±, . . . ,y m ] 

denote the set of SOS polynomials in m variables yi,...,y m , i.e. the set of polynomials that can be 

t 

represented as p = J^Ph Pi e ^m, where P m is the polynomial ring of m variables with real coefficients. 

i=i 

Then, a sufficient condition for (32) is that there exist SOS polynomials r , n, G E [x] and polynomials 
pi, such that 

ro + ■ T ^ + H ■ ■ T vfifi + J2, P lhl + (II 9kf = 
Matlab toolboxes SOSTOOLS [47], or YALMIP [30] automate the process of converting an SOS problem 
to an SDP, which is subsequently solved by available software packages such as LMILAB [18], or SeDumi 
[53]. Interested readers are referred to [42], [35], [43], [47] for more details. 

B. Optimization of Lyapunov Invariants for Mixed-Integer Linear Models 

Natural Lyapunov invariant candidates for MILMs are quadratic and affine functionals. 

I) Quadratic Invariants: The linear parameterization of the space of quadratic functionals mapping 

R n to R is given by: 



Vl = {V : W 1 ^ R | V(x) = 





T 






X 


P 


X 


, P e §™ +1 J 


1 _ 




. 1 _ 





(33) 



where § n is the set of n-by-n symmetric matrices. We have the following lemma. 
Lemma I: Consider a program V and its MILM S (F, H, X , n, n w , n v ) . The program admits a quadratic 

(9, pi) -Lyapunov invariant V G Vj, if there exists a matrix Y G M. neXnH , n e = n + n w + n v + 1, a 

diagonal matrix D v G D n ", a positive semidefinite diagonal matrix D xw G Jf^ nw ) and a symmetric 
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matrix P e S n+1 , satisfying the following LMIs: 

L[PL X - QL T 2 PL 2 < He (YH) + LjD xw L 3 + L*D V L 4 - (A + /i) L?L 5 
X = Trace D xw + Trace D v 

where 















T 


0(n+n„)xn, 


T 




Li = 


F 








, L 4 = 


In v 


; i<5 = 


0(„ e -l)xl 




L 5 




Olx(n e -l) 1 




0(n„+l)x(ra+n TO ) 




_0lXTl„ 




1 



Proof: The proof is presented in Appendix II ■ 

The following theorem summarizes our results for verification of absence of overflow and/or FTT for 

MILMs. The result follows from Lemma 1 and Corollary 1 with q — 2, h — f, though the theorem is 

presented without a detailed proof. 

Theorem 2: Optimization-Based MILM Verification. Let a : -< a ^ I n be a diagonal positive 

definite matrix specifying the overflow limit. An overflow runtime error does not occur during any 

execution of V if there exist matrices Yj € f neXnff , diagonal matrices D iv e B) nv , positive semidefinite 

diagonal matrices D ixw e D" +ril ", and a symmetric matrix P e § n+1 satisfying the following LMIs: 

U l] p [^ l] T <0, VxoGXo (34) 
L\PL\ - 6L T 2 PL 2 ■< He (y^) + L^D lxw L 3 + ^^1^4 - (Ai + fi) LjL 5 (35) 
LfALi - L T 2 PL 2 ■< He (y 2 //) + LjD 2xw L 3 + L T A D 2v U - A 2 ^L 5 (36) 

where A = diag {a~ 2 , — 1} , A, = Trace + Trace D iv , and ^ A™, i = 1,2. In addition, if 
H + 9 > 1, then P terminates in a most T M steps where T u is given in (22). 

2) Affine Invariants: Affine Lyapunov invariants can often establish strong properties, e.g., bound- 
edness, for variables with simple uncoupled dynamics (e.g. counters) at a low computational cost. For 
variables with more complicated dynamics, affine invariants may simply establish sign-invariance (e.g., 
Xi > 0) or more generally, upper or lower bounds on some linear combination of certain variables. As 
we will observe in Section V, establishing these simple behavioral properties is important as they can 
be recursively added to the model (e.g., the matrix H in a MILM, or the invariant sets Xi in a graph 
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model) to improve the chances of success in proving stronger properties via higher order invariants. The 
linear parameterization of the subspace of linear functionals mapping IR n to R, is given by: 

V x = [y : W 1 ^ R | V{x) = K T [x if , K e R n+1 } . (37) 
It is possible to search for the affine invariants via semidefinite programming or linear programming. 

Proposition 5: SDP Characterization of Linear Invariants: There exists a (9, yu)-Lyapunov invariant 
V E Vl for a program V = S (F, H, X ,n,n w ,n v ) , if there exists a matrix F G ]R rieXn » ) a diagonal 
matrix D v 6D n ',a positive semidefinite diagonal matrix D xw G B { + +nw)x{n+nw) , an d a matrix K G R n+1 
satisfying the following LMI: 

Re(LjKL 5 - 9LlK T L 2 ) ■< He(Fif) + LjD xw L 3 + L^D V L A - (A + /i) L^L 5 (38) 
where A = Trace D xw + Trace D v and ^ D xw . 

Proposition 6: LP Characterization of Linear Invariants: There exists a (9, /z)-Lyapunov invariant 
for a program "P = S (F, H, X ,n,n w ,n v ) in the class V*, if there exists a matrix F G IR lxnff , and 
nonnegative matrices G IR lxn ", Z},^, G ]R lx ( n+n ™) ) and a matrix K G IR n+1 satisfying: 

K T L 1 - 9K T L 2 -YH- (D xw - D XW )L 3 - (D v - D V )L 4 - + //) L 5 = (39a) 

£>! + (D v + l r + (D xw + l n+n ^ < (39b) 

D v , D v , D xw , D xw > (39c) 
where D x is either or —1. As a special case of (39), a subset of all the affine invariants is characterized 
by the set of all solutions of the following system of linear equations: 

K T Li - 9K T L 2 + L 5 = (40) 

Remark 1: When the objective is to establish properties of the form Kx > a for a fixed K, (e.g., when 
establishing sign-invariance for certain variables), matrix K in (38)— (40) is fixed and thus one can make 
9 a decision variable subject to 9 > 0. Exploiting this convexity is extremely helpful for successfully 
establishing such properties. 



23 

The advantage of using semidefinite programming is that efficient SDP relaxations for treatment of 
binary variables exists, though the computational cost is typically higher than the LP-based approach. 
In contrast, linear programming relaxations of the binary constraints are more involved than the corre- 
sponding SDP relaxations. Two extreme remedies can be readily considered. The first is to relax the 
binary constraints and treat the variables as continuous variables Vi G [— 1, 1] . The second is to consider 
each of the 2 riv different possibilities (one for each vertex of { — 1,1}™") separately. This approach can be 
useful if n v is small, and is otherwise impractical. More sophisticated schemes can be developed based 
on hierarchical linear programming relaxations of binary integer programs [52]. 

C. Optimization of Lyapunov Invariants for Graph Models 

A linear parameterization of the subspace of polynomial functionals with total degree less than or equal 

to d is given by: 

V d x = |V : W 1 ^ R I V(x) = K T Z (x), K e R N , N = ^ ^ j (41) 

where Z (x) is a vector of length ( n ^ d ) , consisting of all monomials of degree less than or equal to d in n 
variables x\, x n . A linear parametrization of Lyapunov invariants for graph models is defined according 
to (10), where for every i G Af, we have cr, (•) G Vt^\ where d (i) is a selected degree bound for <7j (■) . 
Depending on the dynamics of the model, the degree bounds d (i) , and the convex relaxation technique, 
the corresponding optimization problem will become a linear, semidefinite, or SOS optimization problem. 

1 ) Node-wise Polynomial Invariants: We present generic conditions for verification over graph models 
using SOS programming. Although LMI conditions for verification of linear graph models using quadratic 
invariants and the 5-Procedure for relaxation of non-convex constraints can be formulated, we do not 
present them here due to space limitations. Such formulations are presented in the extended report [49], 
along with executable Matlab code in [57]. The following theorem follows from Corollary 2. 

Theorem 3: Optimization-Based Graph Model Verification. Consider a program V, and its graph 
model G {N,£) ■ Let V : Q n H> R, be given by (10), where a { (•) G Vx (<) . Then, the functions a { (•) , 
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i E N define a Lyapunov invariant for V, if for all (i, j, k) E £ we have: 

-<Tj{T^ (x, w)) + 9%(Ji (x) - 11% E E [x, w] subject to (x, w) E ((X t n IlJJ x [-1, n S£ (42) 
Furthermore, P satisfies the unreachability property w.r.t. the collection of sets Xj_, i E Af \ {0} , if there 
exist £j G (0, oo) , % E M\ {0} , such that 

— 00 (x) G S [x] subject to i G Id (43) 
0i (x) - Si E S [x] subject to x G I, (1 X;_, z G A/"\ {0} (44) 

As discussed in Section IV-A2b, the SOS relaxation techniques can be applied for formulating the search 
problem for functions 0j satisfying (42)-(44) as a convex optimization problem. For instance, if 

( (Xt n n£) x [-1, n s* = {(x, w) | /„ (x, w ) > o, ^ (x, w ) = 0} , 

then, (42) can be formulated as an SOS optimization problem of the following form: 

- a j( T ji ( X , W )) + ( X ) - A*J< ~^2 T pfp -^Tpqfpfq ~^2,Plhl E S [X, , S.t. Tp, T pq E E [x, iw] . 

Software packages such as SOSTOOLS [47] or YALMIP [30] can then be used for formulating the SOS 
optimization problems as semidefinite programs. 

V. Case Study 

In this section we apply the framework to the analysis of Program 4 displayed below. 

/ * EuclideanDivision.c * / 

F0 : int IntegerDivision ( int dd, int dr ) 

Fl : {int q = {0}; int r = {dd}; 

F2 : while (r >= dr) { 
F3: q = q+l; 

F4 : r = r - dr; 

Fm : return r; } 

L0 : int main ( int X, int Y ) { 

LI : int rem = {0}; 

L2 : while (Y > 0) { 
L3 : rem = IntegerDivision (X , Y); 

L4 : X = Y; 

L5 : Y = rem; 

Lm : return X; }} 




Program 4: Euclidean Division and its Graph Model 
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Program 4 takes two positive integers X G [1,M] and Ye [1,M] as the input and returns their greatest 
common divisor by implementing the Euclidean Division algorithm. Note that the Main function in 
Program 4 uses the IntegerDivision program (Program 1). 

A. Global Analysis 

A global model can be constructed by embedding the dynamics of the IntegerDivision program 
within the dynamics of Main. A labeled graph model is shown alongside the text of the program. This 
model has a state space X = J\f x [— M, M] 7 , where J\f is the set of nodes as shown in the graph, and the 
global state x = [X, Y, rem, dd, dr, q, r] is an element of the hypercube [— M, M] 7 . A reduced graph 
model can be obtained by combining the effects of consecutive transitions and relabeling the reduced 
graph model accordingly. While analysis of the full graph model is possible, working with a reduced 
model is computationally advantageous. Furthermore, mapping the properties of the reduced graph model 
to the original model is algorithmic. Interested readers may consult [51] for further elaboration on this 
topic. For the graph model of Program 4, a reduced model can be obtained by first eliminating nodes 
Fn, L 4 , L 5 , L 3 , F , Fi, F 3 , F 4 , and L X) (Figure 1 Left) and composing the transition and passport labels. 
Node L 2 can be eliminated as well to obtain a further reduced model with only three nodes: F 2 , L , L M . 
(Figure 1 Right). This is the model that we will analyze. The passport and transition labels associated 
with the reduced model are as follows: 




Fig. 1. Two reduced models of the graph model of Program 4. 
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T F2F2 : x i— > [X, Y, rem, dd, dr, q+1, r — dr] ^F2F2 : x ^ [Y, r, r, Y, r, 0, Y] 
T L0F2 : x ^ [X, Y, 0, X, Y, 0, X] T F2L * : x ^ [Y, r, r, dd, dr, q, r] 

n| 2F2 : {x | 1 < r < dr - 1} IT F2F2 : {x | r > dr} IT F2Ln : {x | r < dr - 1, r < 0} 

Finally, the invariant sets that can be readily included in the graph model (cf. Section II-B2a) are: 

X L0 = {x | M > X, M > Y, X > 1, Y > 1} , X F2 = {x | dd = X, dr = Y} , X hK = {x | Y < 0} . 

We are interested in generating certificates of termination and absence of overflow. First, by recursively 
searching for linear invariants we are able to establish simple lower bounds on all variables in just two 
rounds (the properties established in the each round are added to the model and the next round of search 
begins). For instance, the property X > 1 is established only after Y > 1 is established. These results, 
which were obtained by applying the first part of Theorem 3 (equations (42)-(43) only) with linear 
functionals are summarized in Table II. 



TABLE II 



Property 


q > 


Y > 1 


dr > 1 


rem > 


dd > 1 


X > 1 


r > 


Proven in Round 


I 


I 


I 


I 


II 


II 


II 


°"F 2 (x) = 


-q 


1 - Y 


1 - dr 


—rem 


1 - dd 


1 - X 


— r 


(^F2F2' ^F2F2) 




(1,0) 


(1,0) 


(1,0) 


(0,0) 


(0,0) 


(0,0) 


(^F2F2' ^F2F2) 


(0,0) 


(0,0) 


(0,0) 


(0,0) 


(0,0) 


(0,0) 


(0,0) 



We then add these properties to the node invariant sets to obtain stronger invariants that certify FTT and 
boundedness of all variables in [— M, M]. By applying Theorem 3 and SOS programming using YALMIP 
[30], the following invariants are found 1 (after post-processing, rounding the coefficients, and reverifying): 

a 1F2 (x) = 0.4 (Y - M) (2 + M - r) a 2F . 2 (x) = (q x Y + r) 2 - M 2 

cr 3p 2 (x) = (q + r) 2 - M 2 a 4F2 (x) = 0.1 (Y - M + 5Y x M + Y 2 - 6M 2 ) 

a 5F2 (x) = Y + r - 2M + Y x M - M 2 cr 6F2 (x) = r x Y + Y - M 2 - M 

The properties proven by these invariants are summarized in the Table III. The specifications that the 

'Different choices of polynomial degrees for the Lyapunov invariant function and the multipliers, as well as different choices for 8, fi 
and different rounding schemes lead to different invariants. Note that rounding is not essential. 
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program terminates and that x e [— M, M] 7 for all initial conditions X, Yg [1, M] , could not be established 
in one shot, at least when trying polynomials of degree d < 4. For instance, oif 2 certifies boundedness 
of all the variables except q, while <t 2 f 2 and <t 3F2 which certify boundedness of all variables including q 
do not certify FTT. Furthermore, boundedness of some of the variables is established in round II, relying 
on boundedness properties proven in round I. Given a (x) < (which is found in round I), second 
round verification can be done by searching for a strictly positive polynomial p (x) and a nonnegative 
polynomial q (x) > satisfying: 

q (x) a (x) - p (x) ( (Tx) \ - M 2 ) > 0, T E {T F2F2 , T F2F2 } (45) 

where the inequality (45) is further subject to boundedness properties established in round I, as well as 
the usual passport conditions and basic invariant set conditions. 

TABLE III 



Invariant o-p 2 (x) = 


°"1F 2 ( x ) 


°"2F 2 (x) , o- 3P2 (x) 


°"4F 2 (x) 


°"5F 2 (x) , o- 6 p 2 (x) 


(^F2F2' ^F2F2) 


(1.0) 


(1,0) 


(1.0) 


(1,1) 


(^F2F2' ^F2F2) 


(1,0.8) 


(0,0) 


(1,0.7) 


(1,1) 


Round I: x? < M 2 for x t = 


Y, X, r, dr, rem, dd 


q, Y, dr, rem 


Y, X, r, dr, rem, dd 


Y, dr, rem 


Round II: x 2 < M 2 for x; = 




X,r,dd 




X,r,dd 


Certificate for FTT 


NO 


NO 


NO 


YES, T u = 2M 2 



In conclusion, (j 2 f 2 (x) or f7 3F2 (x) in conjunction with <r 5F2 (x) or a 6F , 2 (x) prove finite-time termination 
of the algorithm, as well as boundedness of all variables within [— M,M] for all initial conditions X, Y e 
[1, M] , for any M > 1. The provable bound on the number of iterations certified by a 5F2 (x) and <t 6 f 2 (x) 
is T u = 2M 2 (Corollary 2). If we settle for more conservative specifications, e.g., x e [— kM,kM] 7 for 
all initial conditions X, Y e [1,M] and sufficiently large k, then it is possible to prove the properties in 
one shot. We show this in the next section. 

B. MIL-GH Model 

For comparison, we also constructed the MIL-GH model associated with the reduced graph in Figure 1 . 
The corresponding matrices are omitted for brevity, but details of the model along with executable Matlab 
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verification codes can be found in [57]. The verification theorem used in this analysis is an extension of 
Theorem 2 to analysis of MIL-GHM for specific numerical values of M, though it is certainly possible to 
perform this modeling and analysis exercise for parametric bounded values of M. The analysis using the 
MIL-GHM is in general more conservative than SOS optimization over the graph model presented earlier. 
This can be attributed to the type of relaxations proposed (similar to those used in Lemma 1) for analysis 
of MILMs and MIL-GHMs. The benefits are simplified analysis at a typically much less computational 
cost. The certificate obtained in this way is a single quadratic function (for each numerical value of 
M), establishing a bound 7 (M) satisfying 7 (M) > (X 2 + Y 2 + rem 2 + dd 2 + dr 2 + q 2 + r 2 ) 1/2 . Table 
IV summarizes the results of this analysis which were performed using both Sedumi 1_3 and LMILAB 
solvers. 

TABLE IV 



M 


10 2 


10 3 


10 4 


10 5 


10 6 


Solver: LMILAB [18]: 7 (M) 


5.99M 


6.34M 


6.43M 


6.49M 


7.05M 


Solver: SEDUMI [53]: 7 (M) 


6.00M 


6.34M 


6.44M 


6.49M 


NAN 


(^F2F2> ^F2F2) 


(MO- 3 ) 


(1,10- 3 ) 


(I.IO- 3 ) 


(1,10-3) 


(1,10-3) 


(^F2F2> ^F2F2) 


(1,10- 3 ) 


(1,10- 3 ) 


(1,10-3) 


(I.IO- 3 ) 


(1,10-3) 


Upperbound on iterations 


T u = 2e4 


T u = 8e4 


T u = 8e5 


T u = 1.5e7 


T u = 3e9 



C. Modular Analysis 

The preceding results were obtained by analysis of a global model which was constructed by embedding 
the internal dynamics of the program's functions within the global dynamics of the Main function. In 
contrast, the idea in modular analysis is to model software as the interconnection of the program's 
"building blocks" or "modules", i.e., functions that interact via a set of global variables. The dynamics 
of the functions are then abstracted via Input/Output behavioral models, typically constituting equality 
and/or inequality constraints relating the input and output variables. In our framework, the invariant sets 
of the terminal nodes of the modules (e.g., the set X K associated with the terminal node F K in Program 
4) provide such I/O model. Thus, richer characterization of the invariant sets of the terminal nodes of the 
modules are desirable. Correctness of each sub-module must be established separately, while correctness 
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of the entire program will be established by verifying the unreachability and termination properties w.r.t. 
the global variables, as well as verifying that a terminal global state will be reached in finite-time. 
This way, the program variables that are private to each function are abstracted away from the global 
dynamics. This approach has the potential to greatly simplify the analysis and improve the scalability 
of the proposed framework as analysis of large size computer programs is undertaken. In this section, 
we apply the framework to modular analysis of Program 4. Detailed analysis of the advantages in terms 
of improving scalability, and the limitations in terms of conservatism the analysis is an important and 
interesting direction of future research. 

The first step is to establish correctness of the IntegerDivision module, for which we obtain 

o- 7F2 (dd,dr,q,r) = (q + r) 2 -M 2 

The function a 7F2 is a (1, 0) -invariant proving boundedness of the state variables of IntegerDivision. 
Subject to boundedness, we obtain the function 

<t 8 F2 (dd, dr, q, r) = 2r - llq - 6Z 

which is a (1, l)-invariant proving termination of IntegerDivision. The invariant set of node F K can 
thus be characterized by 

X M = {(dd,dr,q,r) G [0,M] 4 |r<dr-l} 

The next step is construction of a global model. Given X M , the assignment at L3: 

L3 : rem = IntegerDivision (X , Y) 

can be abstracted by 

rem = W, s.t. W G [0, M] , W < Y - 1, 
allowing for construction of a global model with variables X, Y, and rem, and an external state-dependent 
input W G [0, M] , satisfying W < Y — 1. Finally, the last step is analysis of the global model. We obtain 
the function a 9L2 (X, Y, rem) = Y xM-M 2 , which is (1, l)-invariant proving both FTT and boundedness 
of all variables within [M, M] . 
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VI. Concluding Remarks 
We took a systems-theoretic approach to software analysis, and presented a framework based on convex 
optimization of Lyapunov invariants for verification of a range of important specifications for software 
systems, including finite-time termination and absence of run-time errors such as overflow, out-of-bounds 
array indexing, division-by-zero, and user-defined program assertions. The verification problem is reduced 
to solving a numerical optimization problem, which when feasible, results in a certificate for the desired 
specification. The novelty of the framework, and consequently, the main contributions of this paper are in 
the systematic transfer of Lyapunov functions and the associated computational techniques from control 
systems to software analysis. The presented work can be extended in several directions. These include 
understanding the limitations of modular analysis of programs, perturbation analysis of the Lyapunov 
certificates to quantify robustness with respect to round-off errors, extension to systems with software in 
closed loop with hardware, and adaptation of the framework to specific classes of software. 

Appendix I 

Semialgebraic Set-Valued Abstractions of Commonly-Used Nonlinearities: 
- Trigonometric Functions: 

Abstraction of trigonometric functions can be obtained by approximation of the Taylor series expan- 
sion followed by representation of the absolute error by a static bounded uncertainty. For instance, an 
abstraction of the sin (•) function can be constructed as follows: 



Abstraction of sin (x) 


x G [ 2 ' 2 1 


X e [— 7T, 7r] 


sin (x) G {x + aw w e [—1, 1]} 


a = 0.571 


a = 3.142 


sin (x) € {x — + aw \ w € [—1, 1]} 


a = 0.076 


a = 2.027 



Abstraction of cos (•) is similar. It is also possible to obtain piecewise linear abstractions by first 
approximating the function by a piece-wise linear (PWL) function and then representing the absolute 
error by a bounded uncertainty. Section II-B (Proposition 1) establishes universality of representation of 
generic PWL functions via binary and continuous variables and an algorithmic construction can be found 
in [49]. For instance, if x G [0, 7r/2] then a piecewise linear approximation with absolute error less than 
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0.06 can be constructed in the following way: 

S={(x,v,w)\x = 0.2[(l + v)(l + w 2 ) + (l-v)(3 + w 2 )],(w,v) G [-1, l] 2 x {-1,1}} (46a) 
sm (x) G {Tx E | x E G S} , T : x E ^ 0.45 (l + v)x + (l-v) (0.2a; + 0.2) + 0.06wi (46b) 

- The Sign Function (sgn) and the Absolute Value Function (abs): 

The sign function (sgn(x) = ll[o,oo) (x) — ll(_ 00) o) (x)) may appear explicitly or as an interpretation of 
if-then-else blocks in computer programs (see [49] for more details). A particular abstraction of sgn (•) 
is as follows: sgn(x) G {v \ xv > 0, v G {— 1, 1}}. Note that sgn (0) is equal to 1, while the abstraction 
is multi-valued at zero: sgn(0) G { — 1,1}. The absolute value function can be represented (precisely) 
over [—1,1] in the following way: 

abs (x) = {xv | x — 0.5 (v + w) , (w, v) G [-1, 1] x {-1, 1}} 
More on the systematic construction of function abstractions including those related to floating-point, 
fixed-point, or modulo arithmetic can be found in the report [49]. 

Appendix II 

of Proposition 4: Note that (13)— (15) imply that V is negative-definite along the trajectories of S, 
except possibly for V (x (0)) which can be zero when rj = 0. Let X be any solution of S. Since V is 
uniformly bounded on X, we have: 

-||V|| oo <V(x(0)<0, Vx(t)eX, t>l. 

Now, assume that there exists a sequence X = (x(0), x(l), . . . , x(t), . . . ) of elements from X satisfying 
(1), but not reaching a terminal state in finite time. That is, x(t) £ X^, Vi G Z + . Then, it can be 
verified that if t > T u , where T u is given by (16), we must have: V (x (t)) < — WVW^ , which contradicts 
boundedness of V. ■ 

of Theorem 1: Assume that S has a solution X—(x (0) , x (£_) , ...) , where x (0) G X and 
x(t-) G X-. Let 

7/j = inf V (x) 
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First, we claim that 7^ < max{V(x (£_)), V(x (i_ — 1))} . If h — I, we have x (£_) G /tT 1 and 
7h < (£_)). If h = f, we have x (i_ — 1) G /i -1 (X_) and 7^ < V(x (t- - 1)), hence the claim. 
Now, consider the 9 = 1 case. Since V is monotonically decreasing along solutions of S, we must have: 

7^= inf V(ar) < max {V(x (t_)), V(x (t_ - 1))} < V (x (0)) < supV(x) (47) 

rrefe- 1 ^-) xex 

which contradicts (17). Note that if /i > and h = I, then (47) holds as a strict inequality and we can 
replace (17) with its non-strict version. Next, consider case (I) , for which, V need not be monotonic 
along the trajectories. Partition X into two subsets X and X_ such that X = X U X_ and 

V (x) < ViG X , and V (x) > Vx G X 

Now, assume that S has a solution X— (x (0) , x (t_ ),...) , where x (0) G X and x (£_) G X_. Since 
V (x (0)) > and 9 < 1, we have V (x (t)) < V (x (0)) , Vt > 0. Therefore, 

j h = inf y (x) < max{V(x(t_)),V(x(t_ - 1))} < V(x(0)) < supF(x) 

which contradicts (17). Next, assume that S has a solution X— (x(0) , ...,x , ...) , where x(0) G X_ 
and x (t-) G X_. In this case, regardless of the value of 9, we must have V (x (i)) < 0, Vt, implying that 
Ih < 0, and hence, contradicting (18). Note that if h = I and either // > 0, or 9 > 0, then (18) can be 
replaced with its non-strict version. Finally, consider case (II). Due to (19), V is strictly monotonically 
decreasing along the solutions of S. The rest of the argument is similar to the 9 = 1 case. ■ 

of Corollary 1: It follows from (21) and the definition of X_ that: 

V (x) > sup {Ha- 1 /* (x)\\ q - l} > sup {Wa-'h (x)^ - 1} > 0, Vx G X. (48) 
It then follows from (48) and (20) that: 

inf V (x) > > supV(x) 

xeh-i(X-) xeXo 

Hence, the first statement of the Corollary follows from Theorem 1 . The upperbound on the number of 
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iterations follows from Proposition 4 and the fact that snp xeX ^ x _ uXoo y \ V (x)\ < 1. ■ 

of Corollary 2: The unreachability property follows directly from Theorem 1. The finite time 
termination property holds because it follows from (12), (23) and (26) along with Proposition 4, that 
the maximum number of iterations around every simple cycle C is finite. The upperbound on the number 
of iterations is the sum of the maximum number of iterations over every simple cycle. ■ 

of Lemma 1: Define x e = (x,w,v, 1) T , where x £ [—1, 1]™ , w £ [—1, l] nw , v £ {—1, l} Uv . Recall 

rp 

that (x, 1) = L 2 x e , and that for all x e satisfying Hx e = 0, there holds: (x + , 1) = (Fx e , 1) = L\x e . It 
follows from Proposition 3 that (9) holds if: 

x T e L[PL x x e - 6x T e LlPL 2 x e < -//, s.t. Hx e = 0, L 3 x e £ [-1, l] n+nw , L A x e £ {-1, 1}"" . (49) 

Recall from the S-Procedure ((30) and (31)) that the assertion a (y) < 0, £ [—1,1]" holds if there exists 
nonnegative constants Tj > 0, % — 1, n, such that a (y) < J] Tj (yf — 1) = y T ry — Trace (r) , where 
r = diag(rj) >z 0. Similarly, the assertion a (y) < 0,Vy £ {—1,1}" holds if there exists a diagonal 
matrix fi such that a (y) < \xi (yf — 1) = y 7 fiy — Trace (/i) . Applying these relaxations to (49), we 
obtain sufficient conditions for (49) to hold: 

x r e L\PL x x e -Bx T e L T 2 PL 2 x e < x T e (YH + H T Y T ) x e +x 7 LlD xw L 3 x e +x 7 L 7 D v L 4 x e -/2-Tmce(D xw +D v ) 
Together with ^ D XW) the above condition is equivalent to the LMIs in Lemma 1. ■ 
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